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Abstract. The interstellar gas flow in the inner disk of M3 1 is 
modelled using a new, two dimensional, grid based, hydrody- 
namics code. The potential of the stellar bulge is derived from 
its surface brightness profile. The bulge is assumed to be triax- 
ial and rotating in the same plane as the disk in order to explain 
the twisted nature of M3 1 's central isophotes and the non circu- 
lar gas velocities in the inner disk. Results are compared with 
CO observations and the bulge is found to be a fast rotator with 
a S-band mass-to-light ratio, Tg = 6.5 ± 0.8, and a ratio of 
co-rotation radius to bulge semi-major axis, 1Z= 1.2 ± 0.1, im- 
plying that any dark halo must have a low density core in con- 
tradiction to the predictions of CDM. These conclusions would 
be strengthened by further observations confirming the model's 
off axis CO velocity predictions. 
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1. Introduction 

For many years, there has been ample evidence that bulges 
of barred spiral galaxies (SB bulges) are triaxial Kor- 
mendy (1982). A spiral galaxy cannot have an oblate spheroidal 
bulge if the position angles of the major axes of the bulge, bar 
and disk are all different unless the bulge is tipped with respect 
to the plane of the disk. Given that ob servations of edg e on 
galaxies show that bulges are rather flat ( Kormendy 1981 ), it is 
most likely that the shortest axis is perpendicular to the plane 
of the disk. 

Using stellar kinematic data, Kormendy and Illing- 
worth (1982) provide evidence that many SA bulges (in un- 
barred galaxies) are rotationally flattened, and consistent with 
being oblate spheroids. However, there is also photometric and 
kinematic evidence to show that triaxial bulges may be com- 
mon in unbarred galaxies. In a similar way to barred galaxies, 
a misalignment between the apparent major axes of the bulge 
and the disk is an indicator of triaxiality, as is a twisting of the 
inner bulge isophotes with respect to those of the disk. T hese 
effects have been seen in NGC 2784 ([Bertola et al. 1988b and 
in a sample of 10 unbarred galaxies ( Zaritsky and Lo 1986 ). 



Stronger evidence for triaxiality comes from combining 
photometry with kinematic data. In an non-axisymmetric po- 
tential, the shape of the rotation curve will depend on the 
position of the line of sight and the major axis of the non- 
axisymmetric component. A slowly rising rotation curve or one 
in which a bump of extreme velocities is seen near the centre 
are indications of triaxiality ( Gerhard et al. 1989 ). By com- 
paring the theoretical gas velocity field derived from a triaxial 
potential with the observed kinematic data, one can determine 
whether or not the assumption of triaxiality is justified. Tri- 
axiality has been demonstrated for both NGC 4845 (Ge rhard 
et al. 1989) and for our Galaxy (|Gerhard and Vietri 1986b- 



It is important to know whether triaxiality is common in 
SA bulges, since any asymmetries in the underlying potential 
of the inner galaxy have a large effect on the motions of inter- 
stellar gas and could provide a method of transporting gas into 
galaxy cores. Hence, triaxial bulges could provide clues to the 
fueling of central starbursts and active galactic nuclei. Further- 
more, large scale movement of gas within a spiral galaxy could 
lead to changes in the overall morphology. In particular, a build 
up of gas in the galactic centre can result in the creation of a 
bar or its destruction. 

Triaxial bulges can also help to constrain the masses of dark 
halos in galaxies. Kinematic observations contain the best evi- 
dence for the existence of dark halos. Flat rotation curves im- 
ply that the mass distribution in a galaxy must extend well be- 
yond the observed stellar or gaseous distributions. However, 
in an axisymmetric galaxy it is not possible to determine how 
much of the mass of the galaxy resides in the disk and how 
much in the halo. Many rotation curves can be fitted by models 
ranging from zero-mass disks to 'maximum disks' (van Albada 
and Sancisi 1986), that is disks that are as massive as possible 
whilst ensuring that the halo is not hollow. 

The velocity fields of barred galaxies provide extra infor- 
mation that can be used to break the disk-halo degeneracy 
(Weiner et al. 200C) . Large jumps in velocity are an indica- 



tion of triaxiality (either bar or bulge). These velocity jumps 
can help to determine a mass model for the non-axisymmetric 
component. Assuming a spherical halo, this can then be used 
to constrain the mass-to-light ratio of the disk and the mass 
fraction of the halo. 
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There is both photometric and kinematic evidence for tri- 
axiality in M31, the Andromeda galaxy. Lindblad (1956)| noted 
that the apparent major axes of bulge and disk are misaligned 
and that the central isophotes are twisted with respect to the 
isophotes of the outer disk. He attributed this to the presence of 



Table 1. Observed parameters 



a triaxial bulge rotating in the plane of the disk. Stark (1977) 
showed that the 'Lindblad twist' is generated by a one param- 
eter family of triaxial models of the bulge of M31, where the 
parameter, <fi, is the angle between the major axis of the bulge 
and the line of nodes of the disk. 

In addition, many observations of interstellar gas have de- 
tected anomalously high velocities in both Hi atomic gas 



(Brinks and Burton 1984) and CO molecular gas (Loinard 



et al. 1995) in the central 10' of M3 1 . [Stark and Binney ( 1 994) 
showed that these velocities would occur in a simple, an- 
alytic, rotating barred potential by calculating closed, non- 
intersecting, periodic orbits in the plane of the disk and arguing 
that interstellar gas would largely follow these orbits. 



This paper aims to build on the work of Stark (1977) 



and Stark and Binney (1994) in two ways. The potential of 
the bulge is derived from observations, using the method of 
Stark (1977) and the surface brightness profile of Walterbos 



and Kennicutt (1988) whilst the gas flows are modelled using a 
2D, isothermal, hydrodynamics code. 

Sect. U provides details of the construction of the mass dis- 
tribution of the bulge from the observed surface brightness pro- 
file. Sect. U describes the hydrodynamics code that was used 
to determine the response of the gas in the ISM to the back- 
ground potential. Sect. ^| discusses the results and their com- 
parison with observations. Sect. || comments on the constraints 
on the dark halo implied by these results. Sect. || concludes. 



2. Modelling M31 

The Andromeda galaxy has been studied for over a thousand 
years. It is one of the few galaxies which can be seen with 
the naked eye and was mentioned by Abu I-Husain al-Sufi in 
his Book on the Constellations of the Fixed Stars in 964 AD. 
As the closest large spiral to our own Galaxy, Andromeda pro- 
vides us with an ideal laboratory for probing the structure and 
dynamics of spiral galaxies. Due to its close proximity, high 
resolution photometric observations are available of both the 
stellar and ga seous components. These includ e maps in U, V, B 
and R bands ( Walterbos and Kennicutt 1988 ), Hi (Brinks and 
Burton 1984) and CO ( ^oinard et al. 1995^ Its large inclina- 
tion angle (77°) means that although it is difficult to determine 
accurate positional data on its apparent minor axis, more infor- 
mation is available on the kinematics of the galaxy. This and 
other observational data is summarized in Table 111 

The model of M31 consists of two components: a rotating 
bulge and an axisymmetric component combining the effects 
of both disk and halo. 



Distance of M31* 


D 


690 kpc 


Total B-band absorption* 


A B 


1.3 mag 


Inclination angle* 


e 


77° 


Systemic velocity* 




-315 km s" 1 


Bulge effective radius* 




10' 


Bulge surface brightness* 


h 


22.2 mag arcsec 


Axis ratio of inner isophotes* 




1.54 


Bulge to disk isophote angle** 


v 


10° 


Max. bulge isophote semi-major axis 5 


''max 


11'. 4 



At the distance o f M31, 1' = 200 pc on the major ax is 
f iHodge (1992)|, * |\Valterbos and Kennicutt (1988) 



Burstein and Heiles (1982) and [ van Genderen (1973) 
" iLindblad (1956)1 s IStark (1977 '~ 



2.1. The bulge 



As shown in |Stark (1977J , a surface brightness profile in which 
the isophotes are elliptical can be transformed into a one pa- 
rameter family of volume density profiles in which the density 
is constant on similar, ellipsoidal shells. Stark's procedure re- 
lies on the use of a number of observationally determined vari- 
ables and the input of two other parameters: the appropriate 
mass-to-light ratio, T, and the angle, <j>, between the major axis 
of the bulge and the disk's line of nodes, both of which are 
determined by comparison with gas kinematics. 

The B-band surface brightness profile of the bulge of M3 1 
follows the r* law of |de Vaucouleurs(1958)| : 



I{s) =I e exp(-7.67[(a/r e )*-l]), 

where r e is the effective radius, I e is the surface brightness at 
r e and s is the semi-major axis length of the isophote. Table [j] 
gives the observationally determined values for M3 1 . However, 
to recover the true bulge surface brightness profile, this must be 
corrected for absorption both in M3 1 and in our own Galaxy. 
Using the val ues of As = 0-32 mag for our Galaxy in the direc- 
tion of M3 1 ( [Burstein and Heiles 1982|), A b = 0.98 mag for the 



bulge region of M3 1 ( van Genderen 1973 ) and Mqb = 5.48 for 
the absolute B-band solar magnitude ( Allen 1973 ), the value of 
I e = 22.2 mag arcsec -2 in Walterbos and Kennicutt (1988)| is 
transformed into I e = 289 Lq pc~ 2 . 

To recover an ellipsoidal luminosity distribution from this 
brightness profile, an Abel-type integral equation has to be 
solved. The solution is itself an integral, which is numerically 
integrated using Romberg's method and multiplied by a con- 
stant .B-band mass-to-light ratio, T^, to generate a 3D mass 
density profile. 

The inner isophotes, which are twisted with the respect to 
those of the outer disk, are used to define the spatial extent of 
the bulge. Hence, the bulge model is truncated at an elliptical 
radius r max = 11 '.4, the length of the semi-major axis of the 
largest twisted isophote. For the model of minimum \ 2 ' /N ', this 
corresponds to a bulge semi-major axis of 3.5 kpc. 

Since the hydrodynamic calculations are carried out in 2D, 
only the forces in the plane of the disk are relevant. Therefore 



S. Berman: Hydrodynamic simulations of the triaxial bulge of M3 1 



3 



the force due to the bulge is calculated in the plane z = 0. The 
force at a point (x, y) is found by splitting the bulge into 300 
concentric, triaxial shells and summing the contributions of all 



shells interior to that point. From Hunter et al. (1988), the com- 
ponent of the force acting at (x, y) in the x direction due to an 
ellipsoidal shell of semi-axes <Xj, Cj, where Oj > hi > Cj 
and Cj lies perpendicular to the plane of the disk, is 



dF x 



2 7r G pifli) hi Ci dai 



dn 



^/(af + K )(b'f + K )(cf + K ) dx 



where p(dj) is the tabulated density of the shell with semi- 
major axis length a, and n is the largest positive root of the 
equation 



V 



(bi + K) {aj + n) 



1. 



A similar relationship exists for the forces in the y direction. 

The region inside the central shell is modelled as a homo- 
geneous ellipsoid of density, p c . This central region is split into 
three smaller shells and the volume weighted densities of these 
smaller shells are summed to give p c . 

Table || details the effect that changes in cj) have on the other 
physical variables in the model. Larger values of (j> are associ- 
ated with fatter, shorter and more rapidly tumbling bulges. 



Table 2. Effect of changing 
mass of the bulge 



on the size, shape, speed and 






a* 


b* 


c* 


a/6 


"l.O 


"l.2 


nt.4 


M § 


6° 


25.0 


12.7 


4.80 


2.35 


51.2 


44.3 


37.4 


0.67 


12° 


19.2 


10.4 


6.05 


1.86 


62.4 


51.0 


45.8 


1.11 


18° 


16.8 


10.1 


6.45 


1.66 


70.2 


56.2 


49.3 


1.35 


24° 


15.5 


9.75 


6.65 


1.58 


76.3 


60.4 


52.2 


1.52 


30° 


14.6 


9.40 


6.80 


1.55 


81.3 


63.8 


54.7 


1.65 


36° 


13.9 


8.95 


6.90 


1.55 


85.6 


66.8 


56.8 


1.76 


42° 


13.4 


8.45 


7.00 


1.58 


89.6 


69.5 


58.7 


1.85 


48° 


12.9 


7.75 


7.10 


1.66 


93.4 


72.1 


60.5 


1.94 



* semi axis lengths (arcmin) 

* bulge pattern speed (km s _1 kpc -1 ) for 1Z = 1.0, 1.2 or 1.4 
§ bulge mass per unit mass-to-light ratio (1O 9 M0) 



2.2. The axisymmetric component 

A simple interpretation of observations of gas velocities in 



M3 1 ( Brinks and Burton 1984 ) suggests that the circular speed, 
v c (r), along the line of nodes of the disk is equal to a con- 
stant velocity, vq, at radii greater than fq. At radii less than 
ro, the circular speed falls linearly to zero. The axisymmetric 
component has a mass distribution contrived to ensure that the 
model's circular speed behaves thus. To determine vq and ro, a 
two component straig ht line is fitted to the CO observations of 
Loinard et al. (1995) . The data are folded about X = and the 
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Fig. 1. Circular speed for best fit model. The solid and dot- 
dashed lines are along the bulge major axis whilst the dotted 
lines are along the bulge minor axis 

inside the bulge and will therefore be experiencing severe non- 
axisymmetric motion. The best fit occurs for vq = 248 km s _1 
and ro = 30.0'. This value of «o is deprojected by a factor sin 9 
for use in the calculation of the forces due to the axisymmetric 
component. 

Since the bulge forces have already been calculated, the 
forces due to the axisymmetric component can be found by 
considering that at radii greater than ro, v A + v B — v$ sin~ 2 9, 
where va and vb are the contributions to the circular speed 
along the line of nodes from the axisymmetric component and 
the bulge, respectively. At radii less than ro, the forces due 
to the axisymmetric component are assumed to fall linearly to 
zero. 

Fig. |l] shows the circular speed curve for the model and the 
contributions of both the bulge and the axisymmetric compo- 
nents. 

Sect. [| considers whether and how the axisymmetric com- 
ponent can be broken up into seperate contributions from the 
disk and dark halo. 

3. Hydrodynamic calculations 

The ISM of M3 1 is assumed to be an isothermal gas respond- 



four most central points are excluded from the fit as they are fall 



ing to the fixed stellar potential described above. Cowie (1980) 
argued that an isothermal fluid is a reasonable assumption for a 
cold, dense ensemble of clouds. The sound speed, c g , is taken to 
be the cloud dispersion velocity, ~ 10 km s _1 , and the surface 
density, S, is the average density of the surrounding clouds. 

The response of the gas is determined using a parallel, hy- 
drodynamics code, GALAHAD, based on the shock capturing 
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FS2 algorithm (van Albada and Roberts 1981, van Albada et 
al. 1982). The code solves the discrete, isothermal Euler equa- 
tions with gravitational source terms on a 2D Cartesian grid, 
which represents half of a rotationally symmetric galaxy. The 
80 x 160 grid cells cover 10 x 20 kpc so that each grid cell 
is a square, 125 pc on a side. This is approximately the same 
resolution as that of the CO observations described later. The 
bulge is positioned so that the major axis is aligned with the 
symmetry axis of the grid. 

The FS2 algorithm is well suited to problems in galactic 
hydrodynamics. It is a second-order, flux-splitting scheme with 
no added numerical viscosity, which utilises operator splitting 
to include source terms and dimensional splitting to perform 
calculations in 2D. Discontinuities in the flow are handled by 
a limiter in the FS2 algorithm rather than by jump conditions. 
The self gravity of the gas is neglected in all simulations. 

The isothermal Euler equations are essentially the conser- 
vation laws of mass and momentum in the x and y directions 
for compressible, frictionless fluids plus an equation of state 
defining the pressure as p = Ec^. In a frame rotating with 
the bulge at a pattern speed fl p , and with gravitational source 
terms, the equations are 




place at a steady rate across the galaxy. These phenomena are 
modelled by setting 



R 



where a is the gas recycling co-efficent, So is the initial gas 
surface density and uq and vo are the initial x and y velocities. 
The values of these and the other hydrodynamical and compu- 
tational parameters used in the simulations are listed in Table |[ 

Table 3. Hydrodynamical and computational parameters 



Number of grid cells 


I x J 


80 x 160 


Gas recycling parameter* 


a 


0.3 pc 2 M^ 1 yr 
1 M Q pc" 2 


Initial gas density 


So 


Constant mass radius 


TM 


10 kpc 


Sound speed* 


Cg 


lOkms" 1 


Galaxy radius 


ra 


10 kpc 


Courant number 


C 


0.5 



Athanassoula (1992), 1 Cowie (1980) 



<9U 


dF 












dx 


ay 








u = 




)■•-( 









R, 




G = 




S = ( 




2Q p u 



where u and v are the gas velocities in the x and y directions, 
$ x and $ y are the gravitational potential gradients in the x and 
y directions and R is defined below. 

The pattern speed, Q p , fixes a co-rotation radius, r co , at 
which, on the major axis, the gravitational and centrifugal 
forces are equal and, in the rotating frame, gas there can stand 
still. The co-rotation radius occurs at the distance along the 
major axis at which a '|° ff = 0, where the effective potential 

$ e ff = $ — ^Qpr 2 - Following Athanassoula (1992) , we set 
r co = IZa, where a is the length of the b ulge semi-major axis . 



However, unlike the galaxy models of Athanassoula (1992) , 
M3 1 is a real galaxy, modelled with a somewhat arbitrary bulge 
edge. Hence, it cannot be assumed that the value of 1Z will be 



1.2 ± 0.2 as quoted in Athanassoula (1992) 



During every simulation, the gas gradually loses angu- 
lar momentum and falls inwards to the galactic centre. As in 



Athanassoula (1992), GALAHAD includes a routine to com- 
pensate for this by crudely simulating the interaction of stars 
with the ISM: star formation occurs and gas is removed in ar- 
eas of high density, whereas stellar mass loss is assumed to take 



Simulations start with the gas at a uniform surface density 
So = 1 Mq pc~ 2 and on circular orbits at a speed given by 
(3va{t) with (3 chosen such that (3 2 v\(r) — GM(rM)/rM, 
where tm = 10 kpc and Af(rjvf) is the mass of the model at 

The boundary conditions are set using two rows of ghost 
cells, outside the boundary of the grid, on which the dynamical 
variables are held constant at their initial values. 

The bulge is linearly introduced over half a rotation period 
and the model is evolved to a quasi steady state over three rota- 
tion periods, corresponding to about 0.3 Gyr. A comparison of 
a run at three and four rotation periods shows very little change 
in the velocity fields over much of the galaxy. However, in the 
region in which gas orbits around the Lagrangian points, sig- 
nificant discrepancies are found. Oscillations with amplitude 
^10 km s _1 occur in the total gas velocity with a period of 
70 Myr. These oscillations persist until at least ten rotation pe- 
riods. 

A number of tests were performed to ensure that GALA- 
HAD performs as it should. One simple and effective test for 
codes used in spiral galaxy simulations is to model an axisym- 
metric galaxy to determine how much gas infall there is due to 
numerical viscosity. After 2 Gyr, 98% of the grid still had the 
same gas density as at the start and only 8% of the initial gas 
mass had fallen into the central 2% of the grid. The density in 
the central cells rose from 1 M Q pc~ 2 to 40 Mq pc -2 . 

Shock-tube tests were performed on GALAHAD to test the 
coding of the FS2 algorithm in the absence of source terms. A 
shock-tube test involves setting up gas on either side of a thin 
membrane at a constant density and at zero velocity. Model 
results were compared to the analytical solution for isothermal 
gas in two separate scenarios: the ID shock-tube and the 2D 
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Fig. 2. ID Shock-tube test with sound speed, c = 350 m s _1 . The left hand panel shows the initial density conditions, the central 
panel shows the gas density after 0.01 s and the right hand panel shows the gas velocity after 0.01 s. In all panels, the solid line is 
the analytic solution and each diamond represents a cell in the simulated results. 



oblique shock-tube. The results are presented for the ID test 
in Fig. H and show that GALAHAD adequately replicates the 
analytical solutions. 



4. Simulation results 

Overall, 87 simulation runs were performed within the three 
dimensional parameter space, (<fi, Tg, 1Z). The ranges of the 
parameters are 12° < <j) < 36°, 5.5 < T B < 8.5 and 1.0 < 
K < 1.4. 

The model which most closely fits the observations is de- 
tailed in Table Q The method which was used t o de termine this 



model of minimum %r/N is described in Sect. 4.3 



Table 4. The model of minimum x 2 /-/V 



Semi-major axis a 

Bulge mass Mb u i ge 

Bulge mass-to-light ratio Tb 

Pattern speed Q p 

Rotation ratio 7Z 

Phi angle (j> 



3.5 kpc 

8.1 x 1O 9 M 

6.5 

53.7kms _1 kpc -1 
1.2 

15° 



4.1. Gas flows 

In the plane z = of triaxial or barred potentials, periodic orbits 
exist in two main families, x\ and X2- x\ orbits exist at all radii 
and are elongated along the bulge major axis. X2 orbits exist in- 
side the Inner Lindblad Resonance (ILR) and are oriented per- 
pendicular to the major axis. In the model of minimum x 2 /N, 
these effects can be seen in panel (b) of Fig. ||, where gas near 
the outline of the bulge is moving in orbits elongated along the 
major axis. Near the centre of the bulge, the gas flows along 
orbits elongated along the minor axis. 

When gas falls inwards and passes the ILR, it shifts from x\ 
to lower energy X2 orbits. In doing so, it creates the 'spray' re- 
gion seen as the diagonal gas flows just above and to the right 
and below and to the left of the centre in panel (b) of Fig. 9. 



Since this area contains transient gas which is not on any peri- 
odic orbits, it is an area of very low density and corresponds to 
the low density regions in panel (a) of Fig. || 

As the gas travels through the 'spray' region, it slows down 
and climbs out of the bulge's potential well. At the end of the 
'spray' region, the gas reaches its minimum kinetic energy, 
abruptly changes direction and follows the gas it has encoun- 
tered back down into the potential well. This abrupt change in 
velocity is seen in panel (b) of Fig. || as the peak in the negative 
of the divergence of the velocity, a signature of shocked gas. 
This corresponds to the high density lanes of gas in panel (a) 
of Fig. p|, which form s where the shock causes gas to pile up. 
Athanassoula (1992) shows that these shocks are intimately 
linked to the dust lanes which appear in many barred galaxies. 

4.2. The dust map 

The high density gas appearing in Fig. |] panel (a) suggests the 
existence of dust lanes in M3 1 . Using the value for the mean 
f?-band absorption in the bulge region of Ab = 0.98 mag, (van 
Genderen 1973), the gas density can be calibrated and an ab- 
sorption map for the bulge can be determined. 

Taking the column density of interstellar hydrogen to be 

TY(Htot) = N(H.{) + 2Ar(H 2 ) = 1.9 x \0 2b A v m~ 2 mag" 1 



and A B = l.3A v (Binney and Merrifield 1998), it follows 



that the mean gas density in the bulge region should be 

S Htot = HM pc- 2 . 

Assuming that CO is a reasonable tracer of both interstellar 
hydrogen and dust in M31 (Neininger et al. 1998), SH tot can 
be used to calibrate the gas densities in the bulge region and 
determine the corresponding B-band absorption map, which is 
shown in Fig. Q 



4.3. Comparison of models 

Fig. |] shows the velocity field for a selection of four different 
models in the frame rotating with the bulge. Each model has 
<fi = 22° implying that the semi-major axis of the bulge a = 16' 
for every model. The solid ellipse in each panel marks the ex- 
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a. x (arcmin) b. x (arcmin) c. Distance (arcmin) 



Fig. 3. The model of minimum "X 2 /N. Panel (a) shows as a greyscale the logarithmic density for the central 20', suggesting 
the existence of dust lanes. Panel (b) shows the velocity field for the same region, in the frame rotating with the bulge. The 
greyscale is minus the divergence of the velocity, again highlighting the shocks. The solid ellipse marks the bulge outline in both 
panels. Panel (c) shows and f2 ± k/2 along the line of nodes (solid lines) and the minor axis (dotted lines), Ci and n have been 
determined from numerical derivatives of the potential. The dot-dash line marks the bulge pattern speed, at 54 km s _1 kpc -1 . 




Fig. 4. Absorption map for the bulge of M31. Derived from the gas density field of the model of minimum /N rotated and 
inclined to the observer's viewpoint. The contours and greyscale show the B-band absorption expected from dust. Contours are 
at 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 5.0, 10, 20 and 90 fi-band magnitudes. 
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Fig. 5. Velocity fields for four different runs in the frame rotating with the bulge. The models in panels (a) and (b) are slower 
bulges with 1Z = 1.4 whereas those in panels (c) and (d) are faster with 1Z = 1.0. Panels (a) and (c) are lighter bulges with = 
5.5 and those in panels (b) and (d) are heavier with T b = 7.5. The pattern speeds for the four models are il p = 50, 52, 70 and 77 
km s _1 kpc -1 four panels (a), (b), (c) and (d) respectively. In all four models, <p = 22° and the semi-major axis of the bulge, a 
= 16'. The solid ellipse in each panel marks the bulge and the greyscale is minus the divergence of the velocity, highlighting the 
shocks. Every computational cell is shown. 



tent of the bulge and the greyscale is minus the divergence of 
the velocity, highlighting the shocks. 

The top two panels panels (a) and (b) are the slower bulges 
(1Z = 1.4 and 51 p = 50 and 52 km s _1 kpc -1 ) whereas (c) and 
(d) are faster (1Z = 1.0 and Ci p = 70 and 77 km s _1 kpc" 1 ). The 



left hand panels, (a) and (c), describe lighter bulges (T^ = 5.5) 
and the two on the right, panels (b) and (d) describe heavier 
bulges (T B = 7.5). 

A comparison of the top two panels with the bottom two 
shows the effect of a faster bulge on the shock morphology. 
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A faster bulge has a smaller co-rotation radius and radius of 
the ILR at which gas transfers between x% and X2 orbits. Since 
this transfer of gas between families of periodic orbits is the 
primary cause of the shocks seen in Fig. & we see that fast 
bulges have shocks which occur closer to their centres. This is 
the case irrespective of whether the bulge is heavy or light. The 
faster bulges also have pronounced secondary shocks, which 
emerge near the end of the major axis of the bulge. 

The left hand panels (a) and (c) show the lighter bulges. In 
comparison with the heavier bulges of panels (b) and (d), they 
have weaker shocks. The shocks for the lighter bulges are also 
slightly closer to the centre and slightly straighter than for the 
heavier ones. 



4.4. Comparison with observations 

The observational data u sed to judge the succe ss of the simu- 
lations is CO data from Loinard et al. (1995). The molecular 



gas is likely to be confined to the plane of the disk and follow 
the periodic orbits described above. CO is a good tracer of the 
majority of mo lecular gas in a galax y, particularly in the case 
of Andromeda ( Neininger et al. 1998 ). Hence, it is amenable to 
modelling using isothermal hydrodynamics. 

However, although there have been recent CO maps made 
of most of the Andromeda galaxy (Loinard et al. 1999, Dame 
et al. 1993), there is very little CO data available in the cen- 
tral bulge region. The reason for this deficiency in data could 
be th at the CO is very cold an d difficult to detect in the inner 
disk ( Loinard and Allen 1998 ) or simply not there (Melchior 
et al. 2000). 

The data from ^oinard et al. (1995) were taken using the 
30-m IRAM millimeter radio telescope in 1993. 26 positions 
along the apparent major axis of M31 were surveyed and 
CO was detected in 16 of them with a resolution of 35" x 
2.6 km s . All 16 data points are used in the comparison with 
the models. 

Each model was rotated by <f> and inclined by to match 
the viewing position on the sky. A cut was then taken along the 
line of nodes of the disk, involving an interpolation from off- 
axis positions, and a comparison made with the CO data. Fig. |] 
shows a comparison between the output of the best fit model 
and the CO data. 

The model fits the data quite well except for one point at 
X = -21.6", v\ os = -422.4 km s _1 which cannot be fit using the 
model described above. This anomalous point may be due to 
a molecular cloud orbiting above the plane of the disk. In this 
case, its line of sight velocity would be smaller than if it moved 
in the plane, since it would be partly moving in the z direction. 

To make a quantitative statement about how good a fit a par- 
ticular model is and which model gives the best fit to the data, a 
likelihood function is derived, where the errors in the position 
and velocity of the CO data are assumed to be Gaussian. Given 
that there are errors in both X and v\ os , the difference between 
the data and the model is found by comparing each datum with 
the point in the model which is closest to it. This takes place 




Fig. 6. Gas velocity from the model of minimum x 2 /N on the 
line of nodes of M31. Crosses are CO detections. The dashed 
line is data from the model on the line of nodes after transform- 
ing to the observing frame. Diamonds are the closest points to 
each data point on the dashed line in terms of data error bars. 
The coordinate system is that of Baade and Arp (1964). 



within error space, where all distances are divided by the error 
on X(35") and all velocities by the error on v\ os (2.6 km s^ 1 ). 

This procedure gives a value of x 2 /N for each model and 
a simple method of ascertaining the model of minimum x 2 /N . 
For this model, cj> = 15°, T B = 6.5, TZ = 1.2 and X 2 /N = 4.7 or 
2.7 with or without the far flung data point respectively. Given 
the various assumptions that have gone into this model, it is 
unlikely that values much closer to\ 2 /N = 1 could be attained. 
Table |] provides details of the parameters of this model. 

Contourplots of \ 2 /N are shown in Fig. ^for 1Z = 1.0, 1.2 
and 1 .4. In each plot, the value for the bulge mass-to-light ratio 
is constrained to 0.5 either side of the best fit model value at the 
90% confidence level. However, as the value of 1Z rises from 
1 .0 to 1 .4, <f) becomes less well constrained. This effect happens 
concurrently with a narrower spread of pattern speeds. The re- 
sult is that the pattern speed varies by 3 or 4 km s _1 kpc -1 in 
each of the panels in Fig. [7] in the 90% confidence regions. 

Fig. m shows the value of x 2 /N as a function of and 1Z 
at a constant value of 4> = 18°. Although the plot was only con- 
structed at a single value of <j>, it does indicate that the favoured 
value of 1Z lies between 1.1 and 1.3. This shows that the bulge 
of Andromeda is a fast rotator. 

Combining all of the above results implies that, at the 90% 
confidence level, 11° < < 24°, 5.7 < T B < 7.3 and 51 < 
fip < 55 km s _1 kpc -1 . 

4.5. Off axis predictions 

Unfortunately, the derived values of x 2 /N are not very robust. 
Using the above method, it is possible to vary the value of 
X 2 /N by two simply by taking the result of a particular model 
at four rather than at three bar rotations. This is enough to turn 
a model which fits the data very well into one which is a much 
worse fit. Therefore it becomes much more difficult to deter- 
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Fig. 7. x 2 /N plots for 1Z = 1.0, 1.2 and 1 .4 in panels (a), (b) and (c) respectively. The spurious point has been removed from the 
data sets. Diamonds mark the position of each model. Contours show the 39%, 50%, 75% and 90% confidence levels, whilst the 
greyscale represents the bulge pattern speeds, ft p , for each model. The value of the best fit model in each set is x 2 '/N = 2.9, 2.7 
and 3.3 for 1Z = 1.0, 1.2 and 1.4 respectively. 
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Fig. 8. x 2 /N plot for <j> = 18° using the data set without the 
anomalous point. Diamonds mark the position of each model. 
Contours show the 39%, 50%, 75% and 90% confidence levels, 
whilst the greyscale represents the bulge pattern speeds, il p of 
each model. The value of the best fit model is x 2 /N = 3.1 for 
TZ= 1.2 and T B = 6.5. 



mine whether a particular model is a good representation of 
M3 1 or not and how well the above results can be trusted. 

A much stronger case could be made if more off axis CO 
data was available in the inner regions of M3 1 . In order to an- 
ticipate such observations, I have provided in Fig. ^| a set of 
predictions from the model of minimum x 2 1 /N for off axis gas 
velocities, for comparison with observations. All predictions 
are for cuts parallel to the major axis, taken at 1' intervals from 
0' to 4' inclusive. 




Fig. 9. Gas velocities from the model of minimum x 2 '/N '. Cuts 
are taken at 0', 1', 2', 3' and 4' parallel to the line of nodes of 
M31. The coordinate system is that of Baade and Arp (1964). 
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These results were all calculated using a distance to M31 
of 690 kpc for ease of comparison with previous work. How- 
ever, using recent results from Hipparcos, Stanek and Gar- 
navich (1998) deduce that the distance is 784 kpc. Adopting 
this value means that 1' = 230 pc. Lengths, whether observed 
or deprojected and the mass of the bulge, all scale in proportion 
to the distance whereas the central density scales as the inverse 
square of the distance to M3 1 . 

5. The dark halo 

As mentioned above, the axisymmetric component is assumed 
to comprise a disk and any halo which exists in M3 1 . This was 
done to reduce the number of free parameters in the model 
whilst remaining true to the observations. However, it is ex- 
tremely important to determine what fraction by mass of the 
inner region of M3 1 consists of dark matter. The inner region 
of M3 1 is defined here as a sphere inscribed by the semi-major 
axis of the bulge of the minimum x 2 /N model i.e. a sphere of 
radius 3.5 kpc. 

In order to do this, the total circular speed curve of the 
model of minimum x 2 /N, shown in Fig. [I], is taken to repre- 
sent the combined mass distribution of the bulge, disk and halo. 
Hence, the total mass at 3.5 kpc, M T = 2.2 x 10 10 M Q . The 
contribution of the bulge to this has already been deduced to be 
8.1 x 10 9 Mq, however there remains a degeneracy between 
the contributions of the disk and the halo. 

5.7. Predictions of CDM for the dark halo 

One method of breaking this degeneracy is to represent the halo 
by an NFW profile (Navarro et al. 1996) derived from cosmo- 
logical A -body simulations, 



p(r) 



1 



(r/r s )(l + r/r s 



where p s is a characteristic density, r s = r2oo/c is a char- 
acteristic radius and c is a dimensionless concentration pa- 
rameter. r2oo is the radius within which the mean density of 
the halo is 200 p C rit, where p cr it = 3H 2 /8irG is the criti- 
cal density of the universe, H = hHo is Hubble's constant, 
H = 100 km s -1 Mpc -1 and h is taken to be 0.75. 

The circular velocity curve for an NFW profile is given by 



^halo(V) 
^200 



1 ln(l + cx) - (cx)/(l+cx) 

X 



ln(l 



c/(l 



where V200 is the velocity at r2oo an d x = r/ r2oo ■ From Navarro 
et al. (1997), at a redshift z = 0, 



^200 



fef200 

km s _1 



kpc. 



Hence, for a given value of Hubble's constant, the halo cir- 
cular speed curve is fully specified by ^200 and c. Navarro 
et al. (1997) find a relationship between the characteristic den- 
sity and radius (or, equivalently, between c and f 200) of a halo 



in a given cosmology, thereby reducing the NFW profiles to a 
one parameter family. 

If V200 = vq = 257 km s _1 , then the method of Navarro 
et al. (1997) gives a halo mass fraction at 3.5 kpc of 58% 
and a disk mass of only 5.8 x 10 9 M Q . However, Navarro 
et al. (1996) state that the maximum rotation velocity can be 
as high as 1.4w 2 oo- If v 200 = «o/1.4 = 184 km s _1 , the halo 
mass fraction at 3.5 kpc falls to 44% and the disk mass rises 
to 7.8 x 10 9 Mq. But in either case, the inner region of M31 
would be heavily dark matter dominated if NFW halo profiles 
are to be believed. 

5.2. Fast bars mean minimal halos 



Debattista and Sellwood (2000) convincingly demonstrate that 
fast, long lived bars cannot co-exist with massive dark halos. 
Their argument is based on the idea that dynamical friction 
from a dense, dark matter halo slows down any rapidly rotat- 
ing bar. Hence, any old bar which is found to be a fast rotator 
precludes the existence of a dark halo, other than one with the 
lowest possible central density consistent with not being hol- 
low. 

'Fast' bars have 1.0 < 1Z < 1.4 and, of the few galax- 
ies where 1Z has been measured, all have fallen within this 
range implying that all barred galaxies have minimal halos and, 
hence, maximum disks. Further, unless the distribution of halo 
densities is bimodal, this conclusion of minimum halos should 
be extended to all spiral galaxies. 

Since the model of minimum x 2 /A f or M3 1 has TZ=\.2 
and the bulge of M31 is older than 6 Gyr ( Trager et al. 2000 ), 
M3 1 must have a low density dark matter halo and a maximum 
disk. The exact value of the halo mass fraction inside 3.5 kpc 
is dependent on the profile adopted for disk and h alo, even for 
a maximum disk. The maximum disk model of Kent (1989) 



implies a halo mass fraction of just 0.5% at 3.5 kpc and, al- 
though the model of minimum x 2 /N is unable to make such a 
quantitative prediction, it is clearly at odds with the NFW halo 
predicted by CDM. 

6. Conclusions 

The mass distribution for the bulge of M31 has been deter- 
mined from observations of the bulge surface brightness pro- 
file. The bulge is triaxial and rotates about the same axis as the 
disk. The GALAHAD code, based on the FS2 algorithm, has 
been used to reproduce observations of CO position and veloc- 
ity in the bulge of the unbarred Andromeda galaxy, particularly 
the extreme velocities which appear in the central 10'. 

By matching the model gas velocity field to the obser- 
vations, it has been possible to constrain the B-band bulge 
mass-to-light ratio, Tg, and the bulge pattern speed, £l p , at the 
90% confidence level, to 5.7 < T B < 7.3 and 51 < tt p < 
55 kms -1 kpc -1 . The angle between the major axis of the 
bulge and the line of nodes of the disk has been constrained to 
11° < (j) < 24°. For the model of minimum x 2 /A, T B = 6.5, 
TZ = 1.2, the semi-major axis a = 3.5 kpc and <j> = 15°. 
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The procedure described above improves on previous bulge 
models by constructing the mass distribution from observa- 
tions. Further, the use of a hydrodynamics code to determine 
the response of the gas to the underlying potential provides a 
more accurate representation of the gas velocity field than cal- 
culations of periodic orbits. However, the oscillations inherent 
in a hydrodynamic model of a galaxy limit the accuracy of any 
simulation and may lead to more slack in the predictions than 
the formal errors imply. 

If the bulge really is triaxial, it should be possible to con- 
strain its parameters further by attempting to observe the off- 
set dust lanes that should accompany the gas shocks described 
above. An absorption map has been provided to assist dust lane 
observations. 

Since the bu lge of M3 1 has been shown to be a fast rotator, 
it follows from |Debattista and Sellwood (2000)| that M31 has 



a maximum disk and a halo mass fraction within 3.5 kpc of 
just a few percent. Th is is in stark contrast to the predictions of 
|Navarro et al. (1996) , which imply that the halo mass fraction 
within 3.5 kpc is between 44% and 58%. 

These conclusions would be strengthened by the confirma- 
tion of the off axis gas velocity predictions which have been 
made above. This could be done by confronting the model with 
further CO observations of the central 15' of M31. 

Given that triaxiality has now been demonstrated in both 
our Galaxy and our nearest large neighbour, M31, it seems 
likely that triaxial bulges are a common feature of both SA and 
SB spiral galaxies. 
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